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To analyse ultrarelativistic nuclear interactions, usually either dynamical models like the string model are employed, or a 
thermal treatment based on hadrons or quarks is applied. String models encounter problems due to high string densities, 
thermal approaches are too simplistic considering only average distributions, ignoring fluctuations. We propose a completely 
new approach, providing a link between the two treatments, and avoiding their main shortcomings: based on the string model, 
connected regions of high energy density are identified for single events, such regions referred to as quark matter droplets. 
Each individual droplet hadronizes instantaneously according to the available n-body phase space. Due to the huge number of 
possible hadron configurations, special Monte Carlo techniques have been developed to calculate this disintegration. 
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Studying nuclear collisions at ultrarelativistic energies 
(-Ecms/micleon ^> 1 GeV) is motivated mainly by the ex- 
pectation that a thermalized system of quarks and gluons 
(quark-gluon plasma) is created [1]. There are essentially 
two directions for modelling such interactions: dynamical 
and thermal approaches. The former ones refer to string 
models [2-7] or related methods [8], supplemented by 
semihard interactions at very high energies [9-12]. Here, 
a well established treatment of hadron-hadron scattering, 
based on Pomerons and AGK rules [13], is extended to 
nuclear interactions. Thermal methods [14-19] amount 
to assuming thermalization after some initial time To, 
with evolution and hadronization being mostly based on 
ideal gas assumptions. 
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FIG. I. High density regions in the x-£ plane for a 
typical event. 

Both methods have serious theoretical drawbacks. 
Even for nuclei as light as sulfur the string models pro- 



duce particle densities that high that the hadrons are 
overlapping. So the independent string model is certainly 
too simplistic, and also considering secondary interac- 
tions as binary collisions among hadrons can theoretically 
not be justified. On the other hand it is also not real- 
istic to consider a homogeneous plasma occupying the 
whole available volume, what is assumed in thermody- 
namic models. To illustrate this, we show in fig. 1 a 
typical event of a string model simulation. We consider 
a central S+S collision at 200 GeV (E cms /A « 10 GeV), 
the transverse coordinates being x and y, the longitudinal 
one (= beam direction) being z; it is useful to consider 
space-time rapidity £ = 0.5 hi(t + z)/(t — z) rather than z. 
In the figure, the hatched regions represent high energy 
density (e > s c = 1 GeV) in the x — ( plane (y = 0). We 
find a couple of intermediate size regions of high energy 
density, representing rest masses of few GeV up to few 
tens of GeV. This demonstrates that neither the inde- 
pendent string model is correct, since these high density 
regions cannot possibly be treated as independent had- 
rons nor the thermal approaches, since we do not have 
one big high density object but rather a couple of medium 
size objects in addition to plenty of ordinary hadrons and 
resonances in particular in the periphery. 

In this paper we introduce a completely new approach, 
more realistic than the string model and more realistic 
than thermal approaches, providing a link between the 
two. Based on the string model, we first determine con- 
nected regions of high energy density (e > s c , for a given 
e c ). These regions are referred to as quark matter (QM) 
droplets. For such regions, the initially produced hadrons 
serve only as a mean to produce the proper fluctuations 
in the energy density. Presently, a purely longitudinal ex- 
pansion of the QM droplets is assumed. Once the energy 
density falls beyond s c , the droplet D decays instantan- 
eously into an n-hadron configuration K = {/ii^2 • • • h n } 
with a probability proportional to £1, with SI representing 
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the microcanonical partition function of an n-hadron sys- 
tem. Due to the huge configuration space, sophisticated 
methods of statistical physics [21,22] have to be employed 
to solve the problem without further approximations. 

Our hadronization scenario is referred to as "microca- 
nonical hadron gas (MHG) scenario" . It is certainly not 
the only one and probably not the most realistic one. 
However, we start with this scenario for a couple of reas- 
ons: the hadron gas scenario is a benchmark, widely used 
in the literature (in a simplified fashion though); the 
MHG scenario can be solved exactly; for massless had- 
rons even an analytical treatment is possible, providing 
very useful tests for the complicated numerical proced- 
ures. After having gained experience in the techniques 
to solve the MHG scenario, we plan to investigate altern- 
atives as well. So the purpose of this paper is not so 
much to promote this particular scenario, but rather to 
show how a dynamical and a statistical treatment can be 
combined. 

The first stage of our approach is the identification of 
high energy density regions, based on the string model, 
which is already discussed elsewhere [20]. Due to the 
empirically found correlation, y = £, between the av- 
erage rapidity y and space-time rapidity £, a hypersur- 
face % T of constant proper time r may be introduced, in 
the central region simply defined by t 2 — z 2 = r 2 . Ap- 
propriate coordinates on % T are the space-time rapidity 
C = 0.51n(t + z)/(t — z) and the transverse coordinates 
x and y. After having used the string model (VENUS 
5.08) to get complete information of hadron trajectories 
in space and time, we may now, for given r, determine 
energy densities on % T and thus locate high density re- 
gions on % T (with e > s c ), as shown in figure 1 for a 
typical example. 

High density regions are considered as QM droplets, 
presently it is assumed that they expand purely longit- 
udinally. Whenever other clusters or hadrons cross their 
way, the two objects fuse to form a new, more energetic 
cluster. Due to the expansion, the energy density of a 
cluster will at some stage drop below s c , which causes an 
instantaneous decay. 

We employ the "microcanonical hadron gas (MHG) 
scenario" for the hadronization: the probability of a 
droplet D — charcterized by the invariant mass E, the 
volume V , and the flavour Q = (Q u , Q d , . . .) — to decay 
into a hadron configuration K = {hi , . . . , h n } of hadrons 
hi is given as 

prob(D K) ~ fi(A') , (1) 

with Sl(A') being the microcanonical partition function 
of an ideal, relativistic gas of the n hadrons hi, 

il(A') = C V ol Cdeg Cident 4> , (2) 

with 

Cmi = (2tt) 3 " ' Cdee = n gi Cident = n ~\ ( 3 ) 



Here, Cd eg accounts for degeneracies (gi is the degeneracy 
of particle i), and Cident accounts for the occurence of 
identical particles in K (n a is the number of particles of 
species a). The last factor 

/n 
\\d 3 Pi 6(E -^6i)6(^pi)6 Q ^ qt (4) 

8 = 1 

is the so-called phase space integral, with Si = \/ m 2 + p 2 
being the energy and fi the 3-momentum of particle i. 
The vector qi = (q™ , qf , . . .) represents the flavour content 
of hadron i. The expression eq. (4) is valid for the centre- 
of-mass frame of the droplet D. 

We have to define a set S of hadron spe- 
cies; we include the pseudoscalar and vector mesons 
(it, K, t], rj' , p, K* , ui, <f>) and the lowest spin-^ and spin- 
| baryons (N, A, S, S, A, £*, S*, SI) and the correspond- 
ing antibaryons. A configuration is then an arbitrary set 
{hi, ...,/}„} with hi G S. 

We are interested in droplet masses from few GeV up 
to 10 3 GeV, corresponding to particle numbers n = \K\ 
between 2 and 10 3 . So we have to deal with a huge 
configuration space, which requires to employ Monte 
Carlo techniques, well known in statistical physics. The 
method at hand is to construct a Markov process, spe- 
cified by an initial configuration A'o, and a transition 
probability matrix p(Kt —> Kt+i)- In generating a se- 
quence A'o, A'i, A'2, • • ., two fundamental issues have to 
be payed attention at: 

• initial transient: starting usually off equilibrium, it 
takes a number of iterations, 7 eq , before one reaches 
equilibrium; 

• autocorrelation in equilibrium: even in equilibrium, 
subsequent configurations, K a and K a +i, are cor- 
related for some range /auto of i. 

In general, both 7 eq and / au to should be as small as pos- 
sible. 

We are going to proceed as follows: for a given droplet 
D with mass E , volume V , and flavour Q, we start from 
some initial configuration A'o, and generate a sequence 
A'o, A'i, . . . , Ki , with 7 eq being sufficiently large to have 
reached equilibrium. If we repeat this procedure many 
times, getting configurations A'j 1 - 1 , Kj 2 \ . . . , these config- 
urations are distributed as Sl(A'). So for our problem, we 
have only to deal with the initial transient, not with the 
autocorrelation in equilibrium. We have to find a trans- 
ition probability p such that it leads to an equilibrium 
distribution Sl(A'), with the initial transient 7 eq being as 
small as possible. 

Sufficient for the appropriate convergence to Sl(A') is 
the detailed balance condition, 

«(A' a ) p(K a -+ Aj) = fi(A'j) p(K b -+ K a ) , (5) 

and ergodicity, which means that for any A' a ,A'j there 
must exist some r with the probability to get in r steps 
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from K a to A'j being nonzero. Henceforth, we use the 
abbreviations 

Cl a := tt(K a ); pab := p(K a —> K b ). (6) 

Following Metropolis [21], we make the ansatz 

Pab = W ah U a b , (7) 

with a so-called proposal matrix w and an acceptance 
matrix u. Detailed balance now reads 
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Uba 
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(8) 



which is obviously fulfilled for 



Uab = F 



&b W ba 
&a W a b 



with some function F fulfilling F(z) j F(z 1 ) = z. 
lowing Metropolis [21], we take 

F(z) = min(z, 1) . 
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(10) 



The power of the method is due to the fact that an arbit- 
rary w may be chosen, in connection with u being given 
by eq. (9). So our task is twofold: we have to develop an 
efficient algorithm to calculate, for given K, the weight 
Sl(A'), and we have to find an appropriate proposal mat- 
rix w which leads to fast convergence (small I eq ). The 
first task can be solved, a detailed description will be 
published soon [23]. In the following we discuss about 
constructing an appropriate matrix w. 

Most natural, though not necessary, is to consider sym- 
metric proposal matrices, w a b = Wb a , which simplifies the 
acceptance matrix to u a b = A(£lj/£l a ). This is usually 
referred to as Metropolis algorithm. Whereas for spin 
systems, it is obvious how to define a symmetric matrix 
w, this is not so clear in our case. We may take spin 
systems as guidance. A configuration K is per def. a set 
of hadrons {hi, . . . , h n } with the ordering not being rel- 
evant, so {7T°, TT°,p} is the same as {p, it , it }. We intro- 
duce "microconfigurations" to be sequences {hi, . . . , h n } 
of hadrons, where the ordering does matter. So for a 
given configuration K a = {hi , . . . , h n } there exist several 



microconfigurations K 



{K 3 



(i)' 



}, with 7Tj 



representing a permutation. The weight of a microcon- 
figuration is 



(11) 



with n a being the number of hadrons of type a. Taking 
for example K = {p,ir° ,ir }, there are three microcon- 
figurations {p, 7r°,7r }, {tt°,p, it } and {7r°,7r°,p}, with 
weight fi(A)/3. 

So far we deal with sequences {hi, . . . , h n } of arbit- 
rary length n, to be compared with spin systems with 



fixed lattice size. We therefore introduce "zeros", i.e. we 
supplement the sequences {hi, . . . , h n } by adding L — n 
zeros, as {hi, . . . , h n , 0, . . . , 0}, to obtain sequences of 
fixed length L. The zeros may be inserted at any place, 
not necessarily at the end. Therefore the weight of a 



microconfiguration K a j with 



relative to the 



without, K a j, is one divided by the number of possib- 
ilities to insert L — n zeros, so from eq. (11) we get 



n 



n\ (L — n)\ 



Cl(Ka) . (12) 



We now have the analogy with a spin system: we have a 
one-dimensional lattice of fixed size L, with each lattice 
site containing either a hadron or a zero. Henceforth, 
we use for microconfigurations with zeros the notation 

(9) 

K a j = {hi, . . . , h^) with h{ being a hadron or zero. 

Since from now on we only consider microconfigura- 
tions with zeros (K a j) rather than configurations (K a ), 
we are going to write K a instead of K a j , keeping in mind 
that a represents a double index, and say "configuration" 
rather than "microconfiguration with zeros" . The ad- 
vantage is that we can use the above formulas specifying 
the Metropolis algorithm without changes. 

We are now in a position to define a symmetric pro- 
posal matrix w(K a —> Kb), with K a = {hf, . . . , and 



Kb = {h\,...,h h L ], 



w(K a —> K b ) 



L(L-l) 



- \ ft 6 hth Av(hth^h\h)) 



(13) 



with 



\V(h?h])\- 




if h\h) <E V{Kh a j ) 
else 



(14) 



where V(hfhj) is the set of all pairs (hihj) with the same 
total flavour as the pair (hfhj). The symbol \V\ refers 
to the number of pairs of V. The term {} in eq. (13) 
makes sure that up to one pair all hadrons in K a and 
Kb are the same, the term 2/L(L — 1) is the probab- 
ility to randomly choose some pair of lattice indices i 
and j. So our proposal matrix amounts to randomly 
choosing a pair in K a , and replacing this pair by some 
pair with the same flavour, with all possible replacements 
having the same weight. The proposal matrix is obvi- 
ously symmetric, since v is symmetric (the symmetry of 
v is crucial!). We have now fully defined an algorithm, 
which due to general theorems will converge, but how 
fast, i.e., how large is 7 eq ? Considering particle ratios, 
like n w o /n w + , we find immediately that we have a very 
slow convergence, so 7 eq is too large for the method to 
be of practical importance. This is obvious, since the 
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FIG. 2. Multiplicity versus the number of iterations. 

method is not very democratic: flavourless particles like 
7T° , p° or also zeros are much more frequently proposed 
than all the rest. This shortcoming can be fixed by de- 
fining w such that two pairs are exchanged rather than 
one, the first pair being replaced by a completely arbit- 
rary pair, the second one by some pair to guarantee fla- 
vour conservation. In addition it is necessary to weight 
the "zeros" differently than the hadrons. This improved 
method violates the symmetry of w, however, the asym- 
metry w a i,/wi, a can be calculated and properly taken into 
account. Further details of the "asymmetric algorithm" 
will be published elsewhere [23]. 

The asymmetric method converges quite fast. As a 
check, we consider massless hadrons, where analytical 
results exist. In fig. 2 we plot the iterated total mul- 
tiplicity N for a droplet with E = 10 GeV and V = 10 
fm 3 , compared to the average multiplicity N from the 
analytical calculation (dashed line). The "equilibration 
time" J eq is roughly given as 

I eq /Nv 10, (15) 

with the initial configuration being two 7To's. 

A major advantage of our method is the fact that it 
is parameterfree. The final distribution of particles from 
the decay of a QM droplet is given by the phase space 
only. The necessary technical parameters J eq and L are 
presently systematically studied. For all types of clusters 
(different E, V, flavour) they have to be chosen as small 
as possible but large enough to not affect the results. 
Preliminary calculations have shown that the method is 
sufficiently fast (1 min per S+S event on a DEC Alpha) 
to be of use for the investigation of ultrarelativistic heavy 
ion collisions. Another systematic study aims at a com- 
parison of our method with a canonical hadron gas, which 
may be useful to optimize L and the starting configura- 
tion. 

Being parameterfree, and with the hadron gas scenario 
representing a benchmark, it will be very interesting to 
simulate nuclear collisions and compare with data. The 
scope of our method, however, is much larger. Our main 
intention is to introduce a completely new way to de- 



scribe ultrarelativistic nuclear collisions, by linking dy- 
namical and thermal approaches, and to show that such 
an approach is technically feasible. 
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